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ABSTRACT 

A new code for following the evolution and emissions of proto-neutron stars during the first minute of 
their lives is developed and tested. The code is one dimensional, fully implicit, and general relativistic. 
Multi-group, multi-flavor neutrino transport is incorporated that makes use of variable Eddington 
factors obtained from a formal solution of the static general relativistic Boltzmann equation with 
linearized scattering terms. The timescales of neutrino emission and spectral evolution obtained using 
the new code are broadly consistent with previous results. Unlike other recent calculations, however, 
the new code predicts that the neutrino-driven wind will be characterized, at least for part of its 
existence, by a neutron excess. This change, potentially consequential for nucleosynthesis in the wind, 
is due to an improved treatment of the charged-current interactions of electron flavored neutrinos and 
anti-neutrinos with nucleons. A comparison is also made between the results obtained using either 
variable Eddington factors or simple equilibrium flux-limited diffusion. The latter approximation, 
which has been frequently used in previous studies of proto-neutron star cooling, accurately describes 
the total neutrino luminosities (to within 10%) for most of the evolution, until the proto-neutron star 
becomes optically thin. 



1. INTRODUCTION 

A proto-neutron star (PNS) is born after the core of a 
massive star collapses to supra- nuclear densities, experi- 
ences core bounce due to the repulsive portion of the nu- 
clear interaction which launches a shock wave that may 
eventually serve to disrupt the entire star in a supernova, 
and leaves behind a compact remnant. The overlying 
star is ejected an d some portion of t he mass may or may 
not fall back fc.f. Uanka et al.ll2007l) . In reality the mass 
of the PNS may increase with time due to this accre- 
tion, but a frequent assumption that is reasonable for 
low mass progenitors, and one adopted here, is that the 
PNS evolves in isolation after the shock has exited. Be- 
cause of the large release of gravitational binding energy 
(2 — 5 x 10 53 ergs), the PNS is initially hot and extended 
compared to a cold neutron star but large portions of 
the mass are still at supra-nuclear densities. Due to the 
high density and reasonably large temperature of this nu- 
clear material, it is opaque to neutrinos of all flavors. In 
the outer regions of the PNS where the density is lower, 
the material is semi-transparent to neutrinos. This hot 
extended object undergoes Kelvin-Helmholtz cooling by 
emitting n eutrinos of all flavors over a period of up to 
a minute ([Burrows fc Lattimerl [19861. at which time it 
transitions to a phase of optically thin neutrino cooling. 

This qualitative description was confirmed when 
about twenty neut r inos were observed from SN 1987A 
(|Bionta et J.lll987t iHirata et al.lll987l ). There were not 
enough events, however, to determine much detail of the 
coolin g process (jLattimer fc Yahilll 19891 : lLoredo &: Lambl 
2002), though limits were placed on the pr operties of 
weakly interacting particles ([Keil et al.l lT997). If a simi- 
lar core collapse supernova were to occur today, modern 
neutrino detectors would see thousands of events. De- 
tailed modeling of the neutrino emission is needed if we 
are to learn about the central engine of core collapse su- 
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pernovae from a nearby event. 

There are numerous other reasons why understanding 
the properties of late-time supernova neutrinos is impor- 
tant, despite the rarity with which they are directly de- 
tected. An important one is the impact of neutrinos on 
supernova nucleosynthesis. Charged current neutrino in- 
teractions in the wind blown from the surface of PNSs de- 
termine the electron fraction of the ejected material and 
thereby constrain its nucleosynthesis. Current uncertain- 
ties in the relative energies of the electron neutrinos and 
anti-neutrinos are large enough to allow for both neutron- 
rich and proton-rich ejecta, w hich may be favorable for re- 
process dWooslev et al.l fl~994[) and i/p-process n ucleosyn- 
thesis ([Frohlich et al.l 120061 : iPruet et all 120061 ). respec- 
tively. Recent work points to the wind ejecta being pro- 
ton r ich at all times (jHudepohl et aLll2010t iFischer et ali 
2010), but the rea sons for this change are only beginning 
to be understood (|Fischer et al.ll20ilt) . Additionally, the 
average energies of [i and r neutrinos also significantly 
affect the neutrino spallation rates that determine nucle - 
osynthetic yields of the i^-process (jWooslev et al.l Fl990). 
which may be responsible for a number of rare isotopes. 

It is also possible that current neutrino detectors with 
upgrades or next generation neutrino detectors will be 
able to observe the diffuse background of neutrinos pro- 
duced by supernovae over the lifetime of the universe 
(|Horiuchi et al.ll2009| ). Predictions for the diffuse MeV 
scale neutrino background density depend significantly 
on the integrated spe ctrum of neutrinos emitted in core - 
collapse supernovae (jWooslev et al.l [19861 : I Andol 120041) . 
The integrated neutrino emission is dominated by PNS 
evolution, so that accurate modeling of PNSs can con- 
tribute to understanding the diffuse supernova neutrino 
background. 

Finally, the neutrino emission from the "photosphere" 
of PNSs gives the initial conditions for the study of both 
matte r-induced and neu trino-induced neutrino oscilla- 
tions (jDuan et al.l 120061) . The differences between the 
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spectra of various neutrino flavors, especially v e and v^r-, 
can significantly affect the impac t of flavor evolut ion in 
the nearly free streaming regime (IKeil et all 12003ft . The 
rate of PNS cooling also has the potential to put h mits 
on exotic physics, such as axions (IKeil et al.lll9"95ft . the 
prese nce of quark matt er or a Kaon condensate in the 
core ([Pons et al.ll20QlaH bT). as well as possible extensions 
of the standard model using data already in hand from 
SN 1987A. 

Theoretical predictions of post-b ounce neutrinos have 
existe d for more tha n 25 ye ars dBurrows fe Lattimerl 
1981 IMavle et al.l fl987t IKeil fc .Tankal 119951: 



Sumivoshi et al.l 119951 iPons et al\ 11 999; 
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the evolution of PNSs is described by the 
Hclmholtz cooling of the collapsed, shock 
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remnant of a core-collapse supernova, it is fundamen- 
tally a radiation hydrodynamics problem (although the 
regions important for neutrino emission are not very 
dynamic after bounce). Over time, the treatment of 
radiative transfer and neutrino microphysics in simu- 
lations has become increasingly sophisticated, moving 
from the equilibrium flux limited diffusion (EFLD) and 
great ly simplified neutrino physics ([Burrows fc Lattimerl 
I1986T) to full solu tions of the Boltzmann equation 
' Fischer et~al\ 12010ft with more realistic microphysics 
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Here, a new fully implicit code is developed for cal- 
culating the detailed evolution of PNSs in spherically 
symmetric general relativity within a variable Edding- 
ton factor formalism. The structure of the paper is as 
follows: In section [3J the equations of neutrino trans- 
port within t he proje c ted sy mmetric trace-free moment 
formalism of IThornel (|1981h are described, and generic 
neutrino source terms for this formalism are derived. In 
section [3l a method for obtaining closure relations for 
the moment equations via a formal solution of the Boltz- 
mann equation are described. A fully implicit numeri- 
cal implementation of neutrino transport coupled to hy- 
drodynamics/hydrostatics is described in section [4] (with 
code tests described in appendix lA"]) . A fiducial model 
of PNS cooling is detailed in section O These results 
are compared with the results of an EFLD calculation 
of PNS cooling in section 16.11 The implications of these 
new calculations of PNS cooling on the composition of 
the neutrino driven wind are discussed in section 16.21 
In section 16. 3| the properties of the integrated neutrino 
emission are discussed. The convention % = c = G = 1 
is adopted in sections [5] through @] to avoid a plethora of 
factors. In section IB\2l units with h = c = 1 are used. 



2. THE MOMENT APPROACH TO GENERAL 
RELATIVISTIC RADIATIVE TRANSFER 

The equations of radiati ve transfer in cu rved space- 
times were first derived by ILindquist] ([1966ft . which de- 
scribed the evolution of the invariant distribution func- 
tion along geodesies in phase space. The general form of 
the general relativistic Boltzmann (or Lindquist) equa- 
tion in the absence of external forces is 

df{x»,p v {x»)) _ p(df_ 1 df_\ = 

dr P \dx? Pi* dp a ) \dr 

ll) 



where / is the invariant distribution function, pP is the 
neutrino four-momentum (which is constrained to be on 
mass shell), and are the Christoffel symbols. The col- 
lision term on the right hand side describes the destruc- 
tion and production of neutrinos on a particular phase- 
space trajectory by capture processes, pair annihilation, 
scattering, and their inverses. In addition to describing 
the propagation of neutrinos along trajectories in phys- 
ical space, this also encodes the evolution of the energy 
of neutrinos along geodesies of the spacetime. In three 
spatial dimensions, this is a seven dimensional equation 
that needs to be solved for each neutrino species. 

A number of numerical strategie s can be employed 
to so lve the transport problem (Mihalas fc Mihalasl 
1984). Foremost among these are discrete ordinate 
methods, where the Boltzmann equation is directly dis- 
cretized i n momentum space as well as in physical 
space (e.g. lYueh fc Buchler1ll977tlMezzacappa fc Messerl 
1999; ILiebendorfer et al.l 12004). and moment-based ap- 
proaches, where angular integrations of the Bol tzmann 
equation in momentum space are performed (e.g. IThornel 
[19811: IBurrows et al.ll2000l:IRampp fc Jankall2002ft . These 
two approaches give similar results in one-dimensional 
models , at least in the context of core-collapse super- 
novae (jLie bcndorfe r et ail [2005) . An additional tech- 
nique that has only been employed for solving static 
problems in the supernova context, but is perhaps the 
most capable of retaining fidelity to the underlying 
Boltzmann equation, is Monte Carlo neutri no transport 
([Janka fc Hillebrandt|[T989l: IKeil et al.ll2003l) . 

The moment approach results in an infinite hierarchy 
of coupled equations which needs to be truncated at some 
order in practice. Generally, only the zeroth and first or- 
der moment equations are retained and a closure relation 
is assumed between the first two moments and the higher 
order moments that enter the first two moment equa- 
tions. Such sche mes are referred to as varia ble Eddington 
factor methods ([Mihalas fc Mihalasl ["l984ft . When only 
the first two moments are used, the number of equations 
relative to discrete ordinate methods is significantly re- 
duced, easing the computational burden (especially in an 
implicit scheme like the one described below). Of course, 
this gain in computational efficiency is useful only if rea- 
sonable closures can be obtained. The closure relations 
only encode information about the angular distribution 
of neutrinos, so that the approximations involved in solv- 
ing a linearized Boltzmann equation do not severely im- 
pact th e fidelity of numerical calculations to th e true so- 
lution ([Mihalas fc Mihalasl [l984t lEnsmanlll994ft . 

Here a variable Eddington factor approach to radia- 
tive transfer is employed, with the closure relations be- 
ing obtained from a formal solution of the static rela- 
tivistic Boltz mann equation . This appr oach is similar to 
that u sed by IBurrows et al.1 ([2000ft and iRampp fc Jankal 
(2002), except for being fully general relativistic, incor- 
porating both inelastic scattering and pair production 
(in contrast to only the former), using energy integrated 
groups rather than energy "pickets" , and in the specific 
method of finding the closure relations. The formalism 
for this method is described below. 

The moments of the Boltzmann equation also most 
naturally give the various forms of the diffusion approx- 
imation, which has been used in the majority of PNS 
studies ([Burrows fc Lattimerl 119861: IKeil fc Ja nka 1995; 
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iPons et g/.| [l999; Robe rts et aT1l2012D and in a significant 
fractio n of studies of the early core-collapse and bounce 
phases (|Bruennlfl985t IWilson fc Mavli |1993[). The for- 
malism is connected to EFLD in appendix [B] 

2.1. General Relativistic Generalities 

In spherical symmetry, it is simplest to work in a co- 
ordinate system that anticipates a Lagrangian frame for 
the fluid. The metric for such a space-time is given by 
(|Misner fc Shardfl96l 



ds 2 = -e 2 Ut 2 + [ - 



da 2 + r 2 dn 2 



(2) 



where ds is the invariant interval, t is the time measured 
at infinity, r is the areal radius, fi is the solid angle, and 
r and <f> are metric potentials. Coordinate freedom can 
be exploited to choose thi s frame to be the rest-fram e of 
the fluid, which demands ()Lieb cndorfcr et al. 20 01aD 

- = 0) 

da A-Kr 2 nB 

Here, n B is the baryon number density and 



„ / 2m . , 

T=^jl + u 2 - — (4) 

where u and m are defined below. With this choice, 
the orthonormal frame associated with the coordinate 
frame is just the rest frame of the fluid and da is just 
the change in enclosed baryon number with the physical 
volume. Therefore, this formulation is working in the 
Lagrangian frame, as claimed. 

The equations of spherically symmetric general rel- 
ativistic hydrodynamics and the Einstein equat ion are 
recorded for convenience isncr fc Shard Il96l . Most 
of these resu lts are nicely presented and detailed in sim- 
ilar form bv iLieben dorfcr ct al. (2001a]). The time evo- 
lution of the areal radius is given by, 

|-A (5) 
which defines u. The evolution of u is given by 



du 



.de* 



Airr 3 (p + Q) 



(G) 



dr r 2 
where Q is the viscosity and p is the pressure of the fluid. 
This gives the equation of hydrostatic balance when the 
left hand side equ als zero (i.e. the Tolman-Oppen heimer- 
Volkov equation (jOppenheimer fc Volkoffl[l93l L The 
enclosed gravitational mass, m, is defined by 



dm / E 
da \Ub 



H 

n B 



(7) 



where e is the internal energy per baryon, E is the total 
neutrino energy density in the rest frame, and H is the 
net radial energy flux from neutrinos. The constraint 
equation for the metric potential <j> is 



e de* 1 
M da nBe^ 



d{e^p) 



1 



<9(r 3 e^Q) 



0, (8) 

e v Oa nse v da r^nse 9 da 
where a small time dependent term has been neglected. 

The transport equations described in the next section 
are formulated in a congruence corresponding to the four- 
velocity field of the PNS (clearly, this is not a geodesic 
congruence). The behavior of this congruence is best 



described by expanding the covariant derivative of the 
four-velocity as 

e 



Uu 



3 ^ 



(9) 



where Up is the tangent four-vector field of the congru- 
ence, P M „ is the projection tensor (which projects into the 
vector subspace orthogonal to U^) a v = U a U" a is the ac- 
celeration, = Uft is the expansion, a^ v is the shear, 
and ujfii, is the rotation. Using the continuity equation, 
the expansion of the congruence becomes 



9 = -£> f ln(ra B ). 



(10) 



In spherical symmetry, the acceleration four-vector is 
parallel to the radial orthonormal basis vector, so that 
only the scalar acceleration is needed 



dr 



(11) 



In spherical symmetry, the shear is characterized by a 
single component, the scalar shear 

»-2S-?0. (12, 
r 3 

Additionally, such a spherically symmetric congruence 
possesses no rotation, so that lj^ — 0. The quantity 



b = 



r 



(13) 



will als o be required, which is related to the extrinsic cur- 
vature (jThornc 198l]). The orthonormal frame temporal 
and radial derivative operators are 

(14) 



D f = e'^ 



Ot 



and 



Df = ^r 2 n B — =T— . 

da dr 



(15) 



2.2. Variable Eddington Factor Transport Equations 

Here, the evolution equations for the neutrino number 
density, energy density, number flux and energy flux are 
derived from the zeroth and first order moments of the 
relativistic Boltzmann equation. The basic results are 
taken from the spherically symmetric version of the pro- 
j ected symmetric trace- free moment formalism of lThornel 
(|1981D . This formalism reduces to an expansion of the 
neutrino distribution function in terms of Legendre poly- 
nomials in a flat space-time. 

The moments of the distribution function are defined 
in spherical symmetry as 



w 



(27T)* 



B n / dfJ,P n ((J,)f{u},(J,) 



where 



B„ 



n\(2n+ 1) 
(2n + l)H 



(16) 



(17) 



P n are the Legendre polynomials, and u is the neutrino 
energy in the fluids rest frame. The first two moment 
equations in sp herical symme try can be read off from 
equation 5.10 of lThornel (fl98lh 

4„ n 3 



w 



d 



+ -@w° + -aw 2 + w\. 



b^w 1 — tt—IjJ 
duj 



© 



-aw 



(18) 
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and 

w 2 f + (a + 3b)w 2 + w\ + 



1 o 4 o 
-w° f + -aw° 



d_ 

duj 



aw + — 



-aw 



-aw 



#9) 



where s are the neutrino source terms defined in sec- 
tion [231 To close this system, define the Eddington like 
factors 

92 = w 2 /w° (20) 

.93 = w^w 1 (21) 

which both go to zero in the limit /(/x) — fo + At/i, 
which corresponds to the diffusion regime. Note that 
these differ f rom the standard defini tion of the Edding- 
ton factors ([Rampp fc Jankal [2002), which is due to 
how I have chosen to calculate the moments. For free 
streaming radiation in a flat background 92 = 2/3, and 
these equations reduce to the linear wave equation for 
h = h(r ± t) = r 2 wo. A method for approximating these 
Eddington factors is detailed in section [3] 

For problems that are close to being static on the ra- 
diation timescale, it is useful to switch the independent 
variable w, the energy in the fluid rest frame, of w l to the 
energy at infinity, v fc.f. iSchinder fc Blud man 198j|). In 
the case of PNS cooling, the energy at infinity is much 
closer to being a constant of the motion and therefore a 
more natural variable. Additionally, this choice simpli- 
fies the formal solution of the Boltzmann equation. The 
moments of the distribution function are then 

w i = w i (r ) v(u, r,t)) (22) 

where the energy at infinity is defined as v = e^'^uj. 
This means that the replacement 

dw l dw l dv dw l 

dx dx dx dv 
needs to be made for all radial and time derivatives, re- 
sulting in 



dw a /riB W 



o 



dt 
e^d_ 
ub dv 



>>B 

e 



92- 2 a 



1 +92 r] w 



d 
da 

v d(j> dw° 



(4nr 2 e< t> w 1 ) 



n dt dv 



and 



~df 



-0 + cr 



w 



n B e 



da 



+ {l~ 92 ) ( a -^ w "~ 



inr 2 e 4 



e 1 "— (24) 
n B 



+ 92\w 



:&g3 w 



E q = du)w , H q = / dujw , 



9 / UWUy , ±±g 

Q"= dojs", and Q] = dws 1 . (26) 

Here, w 9j l is the lower energy bound of an energy group 
and LOgjj is the upper bound. Integrating over energy at 
infinity within groups gives 



and E g = e -0 / dvw° (27) 



N„ 



and similar expressions for F„, H g , and the source terms. 
The operators / dv jv and J dv can then be applied to 
the "red shifted" equations. The evolution the neutrino 
group number densities are described by 



da 



TLB 



e 



-3 + 92 2 a - e dt 



S° 

e*-». (28) 

TIB 



The last term on the left hand side describes the red or 
blue shifting of neutrinos to other groups via compression 
and time variation of the metric potential <j). If the group 
comprises energies from zero to infinity, the red shifting 
terms drop out and one is lef t with the s t andar d number 
transport equation given in iPons et~al\ (|1999f ) . Apply- 
ing the number operator to equation 1251 and simplifying 
gives 



,9Fg 
dt 



e 



F„ 



3e 3 * da 

e 

3" 



:92 



dt 



w 



-Si (29) 



This includes similar terms to equation 1281 plus a term 
that includes the effects of compression on the total num- 
ber flux. The energy group evolution equations are 



d_ 

dt 



riB 



e 

&■ i — 

3 



92 r 



riB 



da 



1 

n B 



e 3 

3 +52 2 a ~ e m 



_^d<t> dw 1 1 
dt 3V (25) 



47rr 2 e 2li 


''H 9 ) 


vv 










n B 




(30) 



To easily deal with optically thick regions where the 
distribution function may possess a sharp Fermi surface, 
energy integrated groups are used rather than discrete 
energy "pickets" . The group numbers, energies, number 
fluxes, energy fluxes, and source terms in group g are 
defined by 

rUs - u dw n f Ug ' u duj i 



Aside from the addition of a compression term and dif- 
ferent factors of e^, this is identical to equation [25J The 
energy flux group evolution equations are 



dt 



:©■ 



3e 4 * da V 91 re* da V 



— W° F 



-W 



93 2 a 



dt 



VW 



•J I 
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The numerical implementation of the red-shifting terms 
is described in section l4~6l 

Additionally, neutrinos have a back-reaction on the 
matter they are propagating through by exchanging en- 
ergy, lepton number, and momentum with the back- 
ground medium. Assuming that the background pos- 
sesses a thermal state, the first law of thermodynamics 
for the medium can be combined with the sum of equa- 
tions [30] over all groups to find an equation for the con- 
servation of total internal energy 



\ 7i ' * An 




3 11 



•£.< 

Si* 



E, 



</■« 



da 



where the sums are over groups and species. Obviously, 
the neutrino energy source terms have exactly canceled 
with the source terms for the medium. 

In the absence of neutrinos, the electron fraction of 
the background medium is fixed, i.e. e~^Y e = 0. 
When neutrinos are included, interactions of electron fla- 
vored neutrinos exchange lepton number with the back- 
ground, yielding e^Y e = ~Ylg^gl n B- The total lep- 
ton number of the medium is given by Yl = Y e + 
Eg \Ng,v e — N g ,p e ] l n B- Combining the evolution equa- 
tion for Y e with equations [28] gives the lepton number 
evolution equation 




N„ 



TIB 



II B 



-— Airr e' 
da \ 



,J E[F g ,u e -F g , Pe } =0. 



(33) 



This constitutes the full set of evolution equations for the 
state of the medium including non-thermal neutrinos of 
all flavors, when the Eddington factors 92 and gs are 
specified. 

2.3. Neutrino Source Terms 

The collision term in equation [1] describes how neu- 
trinos move from one trajectory to another via scat- 
tering and how they are created and destroyed by the 
underlying medium. For the PNS problem these pro- 
cesses include neutral cur rent scattering of f of elec- 
trons, nucleons, and nuclei (jReddv et al.lll998t ). neutrino 
pair production via nucleon-nucleon bremsstrahlung 
(|Hannestad fc Raffeltl Il998f ) and electron-positron anni- 
hilation ([Bruennl I1985D . and charged current processes 
involving electron and anti-electron flavor neutrinos an d 
neutrons and protons, respectively (jReddv et al.lll998l) . 

The details of these microphysical processes are es- 
chewed by assuming that the differential cross-sections 
for these processes are known and referring the reader 
to the paper s cited above, as well as the review 
iBurrows et al.l (|2006l ). The exact details of the micro- 
physics used in the code will be reported in a future 
publication, although certain aspects are discussed in 
sections [5] and 16.21 Many of the results in this se ction 
are well known (e.g. lBruennlli985l ; iPons et «Llll999f ). and 
are included here for completeness and to make clear 



the details of the exact implementation within the inte- 
grated energy group formalism described above. Explicit 
detailed balancing (independent of the choice of under- 
lying scattering kernels) is emphasized. 
The source function for a particular moment is given 

by 

^ l = j^B l f^dfiP l (fi) 
x (j a (i -/)-£+ 3.0- -/)-■£ + - /) - $4) 

which includes contributions from absorption (1/A a ), 
scattering (1/A S ), pair-annihilation (1/A P ), and their in- 
verses (Jajjs, and jp). The choice of metric and reference 
frame implies that the scattering kernels should be evalu- 
ated in the rest frame of the fluid, sim plifying things com- 
pared to hybrid frame approaches (|Hubenv fc Burrows! 
l27)07h . 

For the absorption part, using the standard detailed 
balance relations gives 

ja(l -f)-T~ = TZ (/e?( W > T 'Meq) ~ fl)) , (35) 

where A*" 1 = [1 + cxp{-(w - fi^/T})}^ 1 and f eq is 
a Fermi-Dirac distribution. The scattering contributions 
are 



/ (2^3 W ' 2 / J ^'^(^''^'MO/C^/w) 

(36) 



3s = 
and 

As — 



(2tt) 3 



-1 JO 



J I dfi / dcj)' R s (uj,uj',fi')(l~f (u/, fi ut)), 



and the pair-annihilation contributions are 

3p = 
and 



(37) 



-1 JO 



Ri n {u>,m',fi'){l-f{u'^ out )) 
(38) 



duj' 



(27T) 



-1 JO 



2tt 



A„=/^W 2 / dfi' I d0'R p out (co,co',fi')f(u',fi out ). 



(39) 

T he out g oing co sine is given by fi ou t = fifi' — 
\/l — fi 2 \Jl — fx' 2 cos(j)'. The Rout functions are related 
to the differential cross-section by 



R(w,u/,fi) = 



(2tt) 2 1 da 



(40) 



u/ 2 Vduj'dfi' 

with no phase space blocking term for the final neutrinos 
in the differential cross-section. The R functions obey 
the detailed balance relations for scattering 

R s {oj, w', fi) = R s {lo', w, fi)e {ul ^' )/T . (41) 

and annihilation 

R p out (u,u;',n) = RU^^'^)e {uJ+ -' )/T = R p '(w , J , fi) . 

(42) 

For use in the moment formalism, R must be expanded 
in terms of the Legendre polynomials as 



R{u,w',ii') = Y j R l {u,u')P l { i i') 



(43) 



1=0 
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The distribution function, /, is also expanded in a similar 
way. In general, this results in integrals of the form 

F Mmn = J\d^f\ dn ' f*" d4 >' Pk(jj/)P l (ji)P m (M) 
xP„(^' - \j\ - - i-i' 2 cos 4>') 



and 



2tt 



_£kn_ T. 
2n+l ±lmn 



( 44 ) 

where I lmn = J_ x dliPi(pL)P m (ii)P n {pL). For a 
more detailed description o f such an expansion, see 
iMezzacappa fc Bruennl (|1993l ). 

Using this expansion, the scattering contribution to the 
source term is given by 



4w 3 
W) 



I, ,'2 



B\ I dw uj 



Rjfl -(u-u')/T _ ^ojj 



{21 + If 



21 + 1 



2n+ 1 



a -(u>-u/)/T , > 



(,45) 



m,n— 

and the pair annihilation contribution is given by 
4uj 3 



s l - 

p (2tt) 4 

x[r p S 01 Kl ' 



B t J duj'J 2 e-^+^' T 

Kh mi 



21 + 1 (2i + iy 



1 °° T 

\ " RP f f ' l mn ( 1 

2 ^ ^ n/m/ "2n+l \ 



,(u>W)/T> 



■(46) 



m,n— 



Clearly, all of the moments are coupled to all of the 
other moments by the source terms in addition to the 
coupling present on the LHS of the moment equations. 
Practically, this series must be truncated at some finite 
order. It is standard to use only the zeroth and first mo- 
ment (|Burrows et al.ll2006|) . This convention is followed, 
but with the caveat that this may not be a good approxi- 
mation for the annihilatio n terms near the free streaming 
regime (jPons et al.lll998l) . 

2.3.1. Zeroth Order Source Function 

Now the three contributions to the source functions for 
the number and energy group equations are considered 
separately. Special attention is given to assuring that 
the chosen forms for the source terms explicitly push the 
neutrinos towards equilibrium, independent of the chosen 
opacity functions. 

The absorption part of the source function is given by 



and 



where 



'9 ~ Ng] , 



B„ 



du 



2lo 3 



(2tt) 2 e (^-Mc q )/T + 1 > 



(47) 
(48) 

(49) 



G„ 



duj 



2lo 2 



(27r) 2 e("-^ q )/T + 1 - 



(50) 



The average over the inverse absorption mean free path 
can be performed in a number of ways. For small enough 
energy intervals for groups, the mean free path for the 
central energy of the group can be taken. It can also be 
assumed that the energy within a group is distributed 
as in a blackbody, but renormalized to the total energy 
within the group. Then the averaged inverse mean free 
path is analogous to the Planck mean opacity, but using 
a Fermi-Dirac distribution instead of a Planck distribu- 
tion. Independent of the averaging procedure chosen, 
this term serves to push the neutrino energy density to- 
wards equilibrium with the background medium. 

To find the scattering contribution to the zero order 
moment equation, it is assumed that the distribution of 
energy within a particular energy group is proportional 
to the blackbody distribution. Using this ansatz in equa- 
tion|3S]and then averaging over group energies gives scat- 
tering term 

C=5>o, ffff '^' (j± l) (D g -N g ) 



D„ 



-n, gg >N g [-^-l)(D g ,-N a ,) 



D„ D„ 



and 



9' 

1,99' { r>j 



-K aa 'E g (^-l)(D g ,-N gf ) 



where 



C 



2w 3 f u v 2uj 2 

, d0J j2^' &ndD ^t ^W?- (53) 



The averaged scattering kernel is defined as 



3>i gg , = (i?^V) e - ( ^- )/T ) 



AE„,AE„ 



(54) 



where the average is taken over the energies of the incom- 
ing and outgoing groups. This has the useful property 
$f = &f which derives from the detailed balance 
criterion given above. This form of the scattering source 
term naturally conserves neutrino number, although it 
does not push the neutrino numbers toward the expected 
distribution for a purely scattering process. Both source 
terms go to zero when neutrinos are in both chemical 
and energy equilibrium with the medium. The energy 
exchange expression does not go to zero when g' = g. 
Although it might be naively assumed that this corre- 
sponds to elastic scattering and therefore not contribute 
to the evolution of the group energy, there are in fact 
contributions from small energy transfer scatterings to 
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this term as well. Although these will conserve neu- 
trino number within the group, they can result in energy 
exchange with the medium. This allows the formalism 
to somewhat naturally deal with energy transfer due to 
scattering off nucleons, which generally exchanges energy 
on a scale that is smaller than the group spacing. The 
weighting function for the average over the groups neces- 
sarily involves some level of approximation. The natural 
incorporation of equilibrium far outweighs the small er- 
ror introduced due to the approximate weighting of the 
scattering kernel. 

Using a similar procedure to the one used for the scat- 
tering source term, the energy integrated pair produc- 
tion/annihilation source term is given by 



5 



W-N g ')(D° g -N g ) 



-K gg 'N g N gl {D° g ,/G g , - 1) (D° g /G g ~ 1) 
•*?,,.,-'--,'v (l - e^' T " 



(55) 



where the over bar denotes the energy density and flux of 
a neutrinos anti-species. This term once again naturally 
goes to zero when thermal equilibrium is reached (i.e. 
when N g — G g and F g =0). Here, 

(56) 



*i gg , = {e-^' T R P M"')) 



AE„.AE„ 



The source term Qp can be obtained from the above 
equation by the replacements F g —> H g , N g —> E g , 
D g — > C g , and G g — > B gi while leaving the g' terms 
unchanged. 

2.3.2. First Order Source Function 

For the first order source function terms, one does not 
need to be as careful about getting forms that explic- 
itly go to zero in equilibrium as all terms end up being 
proportional to the first order distribution function and 
therefore satisfy this constraint automatically. 

The absorption contribution to the first order source 
function is 



Si = —F n 



and 



Ql - ~H g 
The scattering source term in equation 



(57) 



(58) 



i.gg' 



F, 



Ng, 


(Dg 


D g >\ 






G g ~) 




(Dg 


Dg>\ 




\G g 


G g >) 



D„ 



Pi 

G 



a) 



The source term Ql can be obtained from the above 
equation by the replacements F g 



H 3 , Ng 



Dg — > Cg, and G g —>B g , while leaving the g' terms un- 
changed. When scattering is iso-energetic, this reduces 
to 

^s.iso— cn ° 



(27T) 



flg(«)-J2f(w)/3 

= -H g [xl{u)- X \{u)/$\ 
where i? ; s is the iso-energetic scattering kernel. 



(60) 



The first order moment equation pair annihilation 
source term is 



g' 



Ng,[l 



(l _ e /3("WA _ D a r 



F, 







D { ; 



> (61) 



and can be found the same replacement required to 
find SI from Q\. 

3. FORMAL SOLUTION OF THE BOLTZMANN 
EQUATION 

To get the factors 32 and 33, an angle dependent ver- 
sion of the Boltzmann equation needs to be solved. First, 
note that the outer layers of the PNS are in tight radia- 
tive equilibrium throughout the duration of the simula- 
tion. Therefore, all time dependence can be reasonably 
dropped in the equation of radiative transfer if one is only 
interested in the ratios of various moments. Of course, 
such an approximation breaks down in highly dynami- 
cal situations. For such circumstances, a c losure scheme 
like the one described in lRampp fc Ja nka (2003) is more 
appropriate. This time-independent formulation of the 
formal solution, which makes calculation of the Edding- 
ton factors sig nificantly easier , is similar to the approach 
advocated by lEnsmanl (|1994l ). except that it incorpo- 
rates gene ral relativistic affects, such a s the bending of 
geodesies. iSchinder fc Bludmanl (|1989ft describe a simi- 
lar, but time-dependent formulation. 

In a spherically symmetri c static spacetim e, the equa- 
tion of radiative transfer is (lLindquistlll966h 



F 



dr \r dr J dfi 



1 



lfe q (r) 



/(v, M, r)] 



+j s [f](l - f(u, n, r)) - A; 1 [f]f(u, /i, r). (62) 

Here, the neutrino distribution function, /, has been 
written in terms of the energy of the neutrinos at in- 
finity, v, and \i is the cosine of the angle of neutrino 
propagation relative to the radial vector. 

A formal solution to equation [62] can easily be found 
using the method of characteristics. The characteristic 
equations are 



dX 



dr 



dfi 



r(i-^) ( i /r _g) df/ V'/a 



(63) 



where A is the physical path length. The second equality 
is easily integrated to find a relationship between r and 
(x along a geodesic. Any geodesic can be characterized 
by the radius at which /1 = 0. First, define the quantity 



(64) 



where the subscript m denotes the minimum radius of 
propagation. This is just the impact parameter of the 
trajectory. Then, for a given /3 and r, the angle of prop- 
agation along a geodesic is given by 



/3e^ 
r 



(65) 
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The first equality in the characteristic equations can 
be integrated to find the physical path length between 
any two radii for a particular characteristic if F and <j> 
are assumed constant over this distance, giving 



aa « ±r~ 



Jrj - e 2 */? 2 - yV| - e 2 ^ 2 



(66) 



where the plus sign is for r/ > r$ and the minus sign 
otherwise. This form is consistent with the assumption 
of constant metric functions across zones (as is used in 
the actual code), but it can introduce difficulties when 
a trajectory moves from one zone to another near the 
radius of minimum propagation. 

Clearly, equation 1621 is a non- linear integro-differential 
equation due to the functional dependence of the scatter- 
ing terms on the local distribution function. An approx- 
imate solution to the Lindquist equation is desired were 
the solutions along characteristics are decoupled and the 
formal solution can be directly integrated. The sim- 
plest approximation is to just to make the replacement 
/ — ► feq i n the scattering terms. This approximation 
will only be valid at high optical depth and is therefore 
suspect for use in the decoupling region. The next order 
approximation is to use a distribution function inferred 
from our knowledge of E g and H g . Assuming that the 
energy is distributed within a group as within a black- 
body gives 

E„ 



fo{y,r) = /eq(^,r) 



fi(v,r) = 3/eq(i/,r) 



(67) 



Employing the Legendre expansion of the scattering ker- 
nel, integrating over outgoing neutrino angle, and only 
including the elastic scattering contribution gives the 
scattering source and sink terms 

Iuj 2 °° 1 

^ ' 1=0 

«/eqH{x5(«)^+^H^} (68) 



and 



(69) 



Using the last characteristic equation, the solution of 
the linearized Boltzmann equation is 

/(f,j8,r/) = /(i/, A r<)e- T(ri,r/) 

+e-^') jf ' ^jf^ { 3s + feJKW) 

where the optical depth is 

r(n,r f ) = j r —(1/AJ + xS). (71) 

This has the appealing property that there is no coupling 
between different /3s and vs, so the evolution of the dis- 
tribution function along each path in phase space can be 
solved for independently. 



4. NUMERICAL IMPLEMENTATION 

Aside from the equation of state and neutrino opac- 
ities for dense matter, PNS evolution is described 
by the transport equations [Ml 121 1301 EH El and 
I3"3l and the structure equations El El El El an d El 

These describe the evolution of the dependent vari- 
ables y(a,t) = {r,u,m,(f>,nB,T,Y e ,Fg,Ng,Eg,H g }. 
To solve these equations numerically, the vari- 
ables _ {a i+ i/2,r i+1 /2, u i+1/2 ,m i+1 /2, F g>i+1/2 , H gti +i/ 2 } 
are discretized on zone edges while the variables 
{(pi,nB,i,Ti,Y Ci i, N g j, Eg i} are discretized on zone cen- 
ters. The derivatives in the PNS evolution equations are 
then finite differenced, turning them to algebraic equa- 
tions for the above dependent variables. This is the most 
natural choice for discretizing the above equations, be- 
cause the thermodynamic quantities, neutrino number 
density, and neutrino energy density are defined on zone 
centers while the neutrino number and energy fluxes are 
defined across zone edges which results in internal energy 
and lepton number conservation being made explicit in 
the discretized equations. 
The general form of these algebraic equations is then 

+(i -9)y (yU,v?, 2/F+0 + o y (y?+i\v? + \yUi) « 

where n is the current time, at which the dependent vari- 
ables are known, and n+ 1 is the next time step at which 
the dependent variables are desired. Here, T denotes 
the differenced time derivatives and y denotes the rest 
of the terms. I choose to employ a fully implicit method 
for solving these equations, i.e. 8 — 1. This leaves a set 
of non-linear algebraic equations that must be solved to 
find the values of the dependent variables at time step 
n + 1. 

These equations are solved by standard high- 
dimensional Ne wton-Raphson (NR) techniques 
(jPress et al.lfl99l) . This requires calculating derivatives 
of all the functions g with respect to y. These deriva- 
tives are calculated analytically. Due to the number 
of derivatives, such an undertaking is prone to error. 
Therefore, all derivative functions are checked against 
numerical derivatives by automated software before 
they are included in the actual evolution code. The NR 
updates are given by the solution of an N z x (6 + 2N g N s ) 
-by-N z x (6 + 2N g N s ) matrix, where N z is the number of 
radial zones, N g is the number of neutrino energy groups, 
and N s is the number of included neutrino species. This 
can rapidly become quite large for reasonable zoning 
and number of energy groups, and become too slow for 
dense matrix techniques. Luckily, the matrix involved 
is in fact block-diagonal, as each zone is only coupled 
to its neighboring zones, which significantly reduces 
computational time compared to solving a general dense 
matrix. 

Although the equations are formally non-linear, they 
are sufficiently close to linear that NR iteration results 
in good convergence after a small number of iterations. 
It is generally demanded that the average relative de- 
viation of the solution from zero is at least less than 
one part in a thousand. Often, the solution found by 
NR iteration satisfies the equations to close to machine 
precision. This scheme has been implemented using 
object-oriented F0RTRAN2003. The block diagonal matrix 
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equations are solved u sing the software package LAPACK 
(jAnderson et al.lll999fl . 

To save computational time, the equations are solved 
using only the neutrino number equations and approx- 
imating the neutrino energy densities and fluxes using 
Eg ss {uj) g N g . Once this set of equations is satisfied, 
a correction step is taken using the energy groups in- 
stead of the number groups. This approximation does 
not seem to introduce any significant error into the cal- 
culation. It is found that the total neutrino energy loss 
calculated using the approximation E g ss (uj) g Ng differs 
from the actual neutrino energy loss by around one part 
in a thousand when thirty energy groups are used. The 
code conserves lepton number to machine precision be- 
cause lepton number conservation is explicitly enforced 
by equation 1331 Conservation of total energy is not ex- 
plicitly enforced. It is found that the total change in 
rest mass over the simulation agrees with the total neu- 
trino energy lost to within a few percent. A series of test 
problems are performed with the code in Appendix [A"l 

4.1. Equation of State 

To close the transport and structure equations de- 
scribed above, an equation of state is required relating 
the pressure, energy density, and equilibrium neutrino 
chemical potential to ns, T, and Y e . Additionally, accu- 
rate derivatives of these quantities are required for calcu- 
lation of the Jacobian matrix for NR iteration. Calls to 
the equation of state must also be computationally effi- 
cient. To meet these requirements, the equation of state 
is implemented in a tabular form. The Helmholtz free 
energy per baryon, F = e — sT, is tabulated as a func- 
tion of tib, T, and Y e , as well as derivatives with respect 
to these variables up to second order. A bi-quintic inter- 
polation is then used to get values of the free energy and 
its d erivatives between grid points (|Timmes fc Swestvl 
2000). This guarantees that the thermodynamic func- 
tions will be smooth in the independent variables, ther- 
modynamically consistent (Swestyj 1996). and does not 
introduce problems in the calculation of the NR correc- 
tions. 

The differential of the Helmholtz free energy is 
P 



dF~- 



-sdT- 



-drii 



ii 



B 



J2 VidY, 



= —sdT+— r dnB + (^ e +Hp-^n)<iY e . (73) 
n B 

From this, the required thermodynamic quantities can 
be read off: 

2 / dF \ / dF s 

P = n B\ "cTi ) , s 



dni 



and/z^cq = (Pe + M P -/•*«) = 



4.2. 



dT 
dF 

dY_ 



(74) 



Neutrino Opacities 

The group averaged neutrino opacities are calculated 
using a ten point quadrature over each group to find an 
effective Planck mean opacity for the absorption terms 
in each group. The scattering and annihilation kernels 
which couple the groups, $g, ff ', are calculated using a five 
point quadrature over both the incoming and outgoing 



energies. Detailed balance is exploited to halve the num- 
ber of calculations required. The scattering terms are 
not weighted by a local thermal neutrino distribution. 

4.3. Integration of the Formal Solution 

The formal solution to the static Boltzmann equation 
is calculated at the beginning of every time step and 
the Eddington factors enter the moment transport equa- 
tions explicitly. Because time independent transport is 
assumed, no previous knowledge of the distribution func- 
tions is required and a new grid of impact parameters 
can be chosen at any time step, withou t having to worry 
about re-mapping old solutions as in iRampp fc Jankal 
l(200l . 

If all quantities are assumed to be constant across 
zones, the formal solution (equation I70p can easily be 
integrated, giving 

f(u, Hs,r i+1 / 2 ) = f(v, /i s , r l _ 1/2 )e _Ari +A/ +A/i (75) 

for the change in / across zone i. 

The physical path length across the zone, AA,, is given 
by equation [6jJ] and the optical depth across the zone is 



At-j = AAj(l/A* + Xo)- 



(76) 



The additions to the neutrino beam from the medium 
and scattering from other beams are given by 



A/ = U 



1 VK + X S 



(1 



and 



H, 



9 n ~&Ti 



AA, 



dXii(X) 



-At, 



5 a(i/a:+x5) 



(77) 



(78) 



The integral required for A/i cannot be calculated an- 
alytically because /z(A) is a fairly complicated function. 
As this is a subdominant term, an "average" /i can be 
pulled out of the integral (which is allowable if /i does 
not change much across the zone). This gives the ap- 
proximation 



1 



ATi 



Bg 1/A* + 1/A* 



(79) 



Note that fi changes most rapidly when it is close to zero, 
but this term contributes the least in that region so the 
error from this approximation should not be too large. 

Numerically, there is a problem with this formulation 
as it stands. Assume that a trajectory in zone i is close to 
its minimum radius of propagation, /i > 0, and that it is 
close to a zone boundary. It then propagates to the zone 
boundary and is considered to be in zone i + 1. Because 
4> is increasing with radius, 4>i+i > 4>i- The new radius is 
taken to be r^^+i, so that the new angle of propagation 
is 



r L,i+l 



-rpi 



(80) 



e <p i+1 -<fn > anc j Tm r L i+1 , Since, //„ must be real, it 
becomes ill defined. In practice this problem is overcome 
setting [x to zero if it would have been imaginary. 

Starting from the outer boundary of the computational 
grid, these equations are solved along an inward going 
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characteristic, through the radius of minimum propaga- 
tion, and then along the outward going characteristic for 
each tangent ray. The impact parameters of the tangent 
ray grid are chosen to be equally spaced in radius for the 
calculations described in this paper. Once the distribu- 
tion function for a particular energy at infinity has been 
calculated along tangent rays, moments of the distribu- 
tion function at radii are calculated from a weighted 
sum that reduces to the correct limit if the distribution 
function is locally constant in angle. Angular resolution 
is reduced at larger depths in the star. Because the dis- 
tribution is extremely close to isotropy and (72 ~ 93 ~ 0, 
this does not pose a significant problem for PNS evolu- 
tion. 

4.4. Boundary Conditions 

To close the system of transport equations, boundary 
conditions for the surface fluxes H g and F g are required. 
For this boundary condition, the formal solution is used 
to calculate the factors 



_ $l.i_dnnf{r,n,v g ) 
/_ x dnf{r, v g ) ' 
a a N a and H a 



so that inbound = otgN g and inbound = a g E g in the 
final zone. At the inner edge of the computational grid, 
incident fluxes are specified (for PNS evolution, they are 
of course specified to be zero). 

The boundary conditions for the radius, gravitational 
mass, velocity, and pressure are implemented by includ- 
ing a fixed ghost zone at the inner and outer bound- 
aries. The boundary condition for the metric potential 
<p is given by matching to the Schwarzschild vacuum so- 
lution to the Einstein equations at the outer boundary. 
This gives <j> s = log(r s ). 

4.5. Rezoning 

To maintain reasonable spatial resolution, conservative 
post time step re-gridding is employed. Where conser- 
vation laws do not specify the properties of a new zone, 
piecewise linear interpolation is used. This generally re- 
sults in smooth radial dependence of the fluid quanti- 
ties. The i mplementation is si milar to the method used 
in Kepler ([Weaver et al.lll978l ). 

The re-gridding is driven by gradients in the density 
and radius. Generally, the radius is not allowed to vary 
by more than 5% between zones and the density is al- 
lowed to vary by no more than 20%. This generally re- 
sults in approximately 100-150 zones being on the grid. 
The choice of relative density changes places high reso- 
lution in regions where neutrino decoupling is occurring. 

4.6. Red Shifting Terms 

Due to red and blue shifting between groups, equa- 
tions [2SJ [2SJ [301 and [3T] contain the un-integrated mo- 
ments w l . Therefore, an approximation method for 
these moments is required. When integrated over all 
energies, these terms go to zero. Therefore, any cho- 
sen numerical scheme must have terms balancing be- 
tween groups for energy conservation. To move forward, 
something must be assumed about how energy is dis- 
tributed in the groups. The simplest scheme is to as- 
sume that it is uniform. Then within a particular group 



w ' 1 = {Ng,F gi E gi H g }/(LUg tH - w b ,l)- It could also be 
assumed that the internal energy is distributed as a black 
body, which is consistent with the assumption used in the 
source terms. The uniform distribution is chosen due to 
its simplicity. For the H g evolution equation, this results 
in 



-Ug,U 
-Ug,L 



(f 



Ho 



H, 



2Au 



s+i 



ff3§er 



-<j>d±\ ( Hg-! 

dt J \2Auig-i 



2Aw, 



(82) 

Similar expressions result for equations |2"51 |2"91 and 1501 It 

is straight forward to verify that these terms disappear 
when summed over groups. 

5. PROTO-NEUTRON STAR EVOLUTION 

Rather than follow the collapse of a massive stellar 
core through bounce, the calculations here start from a 
separate calculation of the highly dynamic phase of ini- 
tial collapse. For ease of comparison with pre vious work, 
the 1. 6 Mq baryonic mass initial model from iPons et~al\ 
(1999) is employed and the affects of convection are not 
considered. This will correspond to a 1.4 Mq gravita- 
tional mass neutron star after it has cooled and can be 
thought of as repres entative of a standard neutron star 
(jKiziltan et al.ll201Q[ ). The cooling and de-leptonization 
of this object is followed for 55 seconds, which is shortly 
after the time the PNS becomes optically thin. 

5.1. Physical Ingredients 

A relativistic mean field of equation of state con- 
sisting of only neutrons, protons, and electrons is 
assumed. T he GM3 parameter set is used with- 
out hyperons (|Glendenning fc Moszkowskilll99lI ). which 



is what was 
opacities are 



used 

also calculated 



Pons et al\ 



(199 



Neutrino 
relativistic mean 



in the 

field approximation using the formalism of iReddv et al.l 
(1998). The tensor polarization is also included so 
that "weak magnetism" affects are included to all or- 
ders ([Horowitz fe Perez-Gar cfa 2003) . The electron scat- 
tering rates from lYueh fc BucMerT([T977l) are used for 
the inelastic scattering kernels. Nucleon scattering is 
assumed to occur within a single group, although the 
opacities are calculated using the full inelastic differen- 
tial cross-sections. Bremsstrahlun g is implemented us- 
ing th e structure function given in lHannestad fc Raffeltl 
(1998). Rather than include this in the annihilation ker- 
nels, the Bremsstrahlung mean free path has been calcu- 
lated assuming a thermal distribution for the secondary 
neutrinos. Given the uncertainty in the Bremsstrahlung 
rate itself and its large density dependence, this is a rea- 
sonable a pproxima tion. Electron positron pair annihila- 
tion (jBruenn| [l985) is also included. Pure neutrino pro- 
cesses (i.e. v e + v e — > u T + u T ) are not included. This 
set of rates is ful ly consistent with the rate set used in 
IPons et all (|1999( 1. but differs significant ly from the rate 
sets used in recent collap se simulations ([Hiidepohl et all 
2010: lFischer et all 120111) . 

The study here uses 30 logarithmically spaced energy 
groups from 2 MeV to 75 MeV plus one final group ex- 
tending from 75 MeV to 1000 MeV to encompass the 
tail of the thermal distribution. This final group is only 
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Figure 1. The internal structure of the PNS for selected times in the fiducial simulation. The temperature, equilibrium electron neutrino 
chemical potential, proton neutron chemical potential difference, dimensionless entropy per baryon, electron fraction, and baryon density 
are plotted at times s, 1 s, 5 s, 10 s, 20 s, and 30 s, from left to right and top to bottom. The horizontal axes show the encl osed baryon 
numbe r in units of the number of baryons in the sun (Nq = 12.04 X 10 56 ). This figure can be directly compared to figure 9 of lPons et al.\ 
(1999), as it was produced using the same initial model and a very similar nuclear equation of state and neutrino opacity set. 



populated deep in the PNS and it is in tight thermal 
equilibrium due to the extremely short mean free paths 
for such high energy neutrinos. Minimal differences are 
found in the PNS evolution if only 20 groups are em- 
ployed to cover the same energy range. 

The adaptive radial gridding algorithm is set to keep 
approximately 130 zones on the grid and allow for at 
most a 20% change in density across a zone and a 10% 
change in radius across a zone. The boundary pressure 
is set so that the outer edge of the model has a density 
around 2 x 10 9 gcm~ 3 . This is a sufficiently low density 
that all of the neutrinos have decoupled well within the 
outer boundary. 

5.2. Structural Evolution 

Qualitatively, the internal structure of the PNS evo- 
lution follows the standard picture of Kelvin-He l mholt z 
PNS cooling a s descr ibed bvlBurrows fc Lattimerl (1986), 
iKeil fc Jankal (|1995l) . and lPons et all (|1999l h ^here the 
gravitational binding energy of the compact object pro- 
vides energy lost to neutrino emission. After the shock 
produced by the supra- nuclear density bounce of the core 
propagates through the outer layers of the PNS, a high 
entropy shocked region is left on top of a cold un-shocked 
PNS core, which has an entropy similar to the initial en- 
tropy of the pre-supernova iron core. The outer shocked 



layers have de-leptonized during the v e burst, but neu- 
trinos in the core itself have been trapped since before 
bounce (although partial deleptonization has occurred), 
resulting in a large non- zero /^^oq and Y e f» 0.3 (c.f. 
ILiebendorfer et aLll2001bl ). This provides the initial con- 
dition for PNS cooling. 

The internal structure of the PNS simulation is shown 
in figure [1] for a number of times (with time zero corre- 
sponding to the starting point of the simulations, not the 
time of core-bounce). The models start with a core en- 
tropy of ~ 1.2. The entropy rises from 1.6 at an enclosed 
baryonic mass of ~ 0.6 Mq to 7.4 at an enclosed mass of 
1.0 Mq. This implies that the supernova shock was born 
at around 0.6 Mq, which is reasonably consisten t with 
the core-collapse results of iThompson et all (|2003). The 
shocked mantle is at low density relative to the core and 
extends to large radius (material that is at a density of 
10~ 5 fm~ 3 is found at 99 km), mainly due to the thermal 
contribution to the pressure. 

From this initial state, the shock heated mantle rapidly 
contracts over the first second or so of the simulation. 
This contraction is driven by the rapid loss of energy 
and lepton number via neutrinos, which can readily es- 
cape due to the low density of the envelope and long 
interaction mean free paths. The loss of lepton num- 
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ber and thermal energy reduces pressure support in the 
mantle, and the mantle responds by rapidly contracting 
(i.e., rapid relative to the cooling timescale of the core, 
not rapid compared to the dynamical timescale of the 
envelope). By two seconds into the simulation, mate- 
rial at a density of 1CP 5 fm~ 3 is at 17 km. This is fairly 
close to the cold neutron star radius for GM3 (13.5 km). 
The work provided by this contraction is enough to in- 
crease the peak temperature of the mantle from 22 MeV 
to 45 MeV even though the entropy of the mantle has 
decreased from 7 to 4 over this period. 

This period of the PNS evolution is most likely to be 
sensitive to the initial conditions for the simulations, as 
at later times the details of the initial structure should be 
washed out. The envelope of the PNS should also be con- 
vective, which significantly alters the rate of energy and 
lepto n number transport in the PNS (c.f. IRoberts et al.l 
120121) . Additionally, there might be significant accre- 
tion luminosity over this period (although this is ap- 
proximately accounted for by the mantle). Therefore, 
especially given the older provenance of the initial con- 
ditions, the results from this period should be taken as 
only qualitatively correct. 

While the mantle is contracting, v e s and v x s are being 
transported down the positive radial temperature gra- 
dient into the core while the v e s are being transported 
outwards down the large equilibrium chemical potential 
gradient. This results in a net heat flux into the core 
and a net lepton flux out of the core. This has been 
referr ed to as "Joule heating" o f the core in previous 
work (Bu rrows fc Lattimerl 119861 ). Additionally, the in- 
ner regions contract over this period due to the increased 
boundary pressure on the un-shocked core from the cool- 
ing mantle. This contributes to the temperature increase 
in the core in addition to the Joule heating. 

After the initial period of mantle contraction, the den- 
sity structure of the PNS becomes similar to that of a 
cold PNS. Joule heating continues to increase the tem- 
perature of the inner most regions until the central tem- 
perature reaches its peak value of 35 MeV at 18 s in the 
simulation. Then, the temperature of the entire star falls 
with time. The entropy evolution exhibits a similar be- 
havior. Lepton number is lost from the entire PNS core 
over this time and the electron fraction evolves toward 
the expected value for matter in beta-equilibrium with 
no net electron neutrino number. After about 15 sec- 
onds, contraction slows since the PNS is nearly at the 
cold neutron star radius. After this, neutrino emission is 
powered chiefly by the loss of thermal energy from the 
star. 

It is also worth noting that the temperature gradient 
and the /x„ c iCq gradient in the shocked layers of the PNS 
become increasingly shallow from 1 s onwards. Addition- 
ally, as the density of the outer layers rises, the neutron 
proton chemical potential difference jX = /j, n — [i v gets 
larger. The increase in /t and the decrease in /i^.cq bring 
p, close to the electron chemical potential \i e = fi + \i Vt iCq 
as time goes on. These considerations have significant 
consequences for the spectral evolution of the neutrinos 
and which are discussed in section 16.21 
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5.3. Emergent Luminosity and Spectral Evolution 

The total integrated energy loss in neutrinos over the 
duration of the simulation is E v = 2.32 x 10 53 erg, and 



Time (s) 

Figure 2. Top panel: Number luminosities as a function of time 
for f e (solid black line), v e (dashed black line), v x (dot-dashed 
black line) and the de-leptonization rate, N Ve —Np e . Bottom panel: 
Energy luminosities as a function of time. The black lines are the 
same as in the top panel, but the solid red line is the total energy 
emitted in neutrinos per time. 



the total lepton number radiated is Nl — 3.2 x 10 56 . 
The neutrino emission from the PNS is shown in figure 
[2j As is discussed above, the first couple of seconds are 
dominated by the contraction of the PNS mantle. Over 
the first second of the simulation, 38% of the total neu- 
trino energy loss and 20% of the total lepton number loss 
occurs. During this period, the v x number luminosity is 
produced mainly by the un-shocked core, as the /x and 
t neutrinos are mainly coupled to the envelope through 
scattering. Therefore, the luminosity in these flavors is 
lower because of the smaller emitting surface (which is 
not offset by the temperature of the core). 

During mantle contraction, there is a high de- 
lcptonization rate driven by the outermost layers of 
the star. After the first few hundred milliseconds, 
de-leptonization slows as the outer layers go towards 
Mfe^q ~ and de-leptonization is driven by diffusion 
out of the core. The values of the de-leptonization rate 
before 80 ms are unrealistic, as they are determined by 
the relaxation of the assumed initial conditions for the 
neutrinos. 

Over the first two seconds, the v x luminosities are sig- 
nificantly lower than the luminosities of the electron fla- 
vored neutrinos. The neutrino energy and number lu- 
minosities as a function of radius at 500 ms after the 
beginning of the simulation are shown in figure |3] First, 
this illustrates that the /j, and r neutrino number fluxes 
are being set much further inside the star (at around 18 
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Figure 3. Neutrino number and energy luminosities at infinity 
as a function of radius at 500 ms into the simulation. The solid 
lines are for electron neutrinos, the dashed lines are for electron 
antineutrinos, and the dot-dashed lines are for [i and r neutrinos. 



km) than the electron neutrinos, but they exchange en- 
ergy out to a significantly larger radius via scattering. 
Second, there is an inward directed anti-electron, /it and 
t flux near the mantle core boundary. As cooling pre- 
cedes, heat diffuses down the positive temperature gra- 
dient (and positive equilibrium chemical potential gradi- 
ent for the anti-electron neutrinos) into the lower entropy 
core. This is the Joule heating discussed above. In con- 
trast, the large negative equilibrium chemical potential 
gradient for the electron neutrinos overwhelms the pos- 
itive radial temperature gradient and the electron neu- 
trino flux is positive everywhere. 

After the PNS has contracted to close to the cold neu- 
tron star radius, the v x luminosity has increased relative 
to the v e and v e luminosities. In fact, the v x luminos- 
ity is about twice the luminosity in either of the electron 
neutrino species. These neutrinos decouple further inside 
the PNS and are therefore emitted at a higher effective 
temperature, resulting in a larger number and energy lu- 
minosity. Between thirty and forty seconds the PNS be- 
comes transparent to neutrinos and the luminosity drops 
off significantly. 

The average energies of the emitted neutrinos at infin- 
ity as a function of time are shown in figure|3J Within the 
integrated energy group formalism, the neutrino energy 
moments at infinity are defined as 

where (uj) g is a group averaged energy and 4> s is the sur- 
face value of the metric potential. 

During the mantle contraction phase, there is the 
standard hierarchy of neutrino average energies (tv c ) < 
(e Pe ) < (e Vx ). After mantle contraction has ceased, the 
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Figure 4. Average energies of the emitted neutrinos measured at 
infinity. 
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Figure 5. Neutrino spectra for all three flavors tracked in the 
simulation. The top panel is at 100 ms after the start of the simu- 
lation, the bottom panel is the spectrum at 5 seconds. 



energy decoupling radius of electron neutrinos and n and 
t neutrinos becomes similar and for the rest of the PNS 
evolution (e Pe ) « (e^)- This is in contrast to the dif- 
ference between the electron neutrino and anti-neutrino 
average energies, which obey {e Va ) < (e Pc ) for the entire 
calculation, although the two average energies get closer 
at late times. An analysis of why this is, its implica- 
tions, and a comparison to other results in the literature 
is given in section f6. 21 

The emitted neutrino spectra at two representative 
times are shown in [5] for reference. The v x neutrinos 
decouple further in the star than the v e neutrinos at five 
seconds into the simulation so that they have a larger 
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number luminosity due to the larger temperatures found 
there, but these flavors have a similar energy sphere 
due to inelastic scattering which accounts for the sim- 
ilar value of the peak of the luminosity as a function of 
neutrino energy. 

6. DISCUSSION 

6.1. Comparison to EFLD 

Un til recently (Tlludcpo hl et al\ l2010t iFischer et all 
2010), most studies of PNS cooling used the 
EFLD approximation to descr i be neutrino trans- 
port (iBurrows fc Lattimerl 119861: iKeil fc Jankal 119951 : 
iPons et aZJI1999URoberts et alJl2012f T It is thus a "worth- 
while exercise to compare the results obtained using 
EFLD and the present variable Eddington factor method 
for transport. No detailed comparison of the effect of 
different flux limiters is attempted since EFLD clearly 
breaks down in the decoupli ng regime ind e pende nt of 
the flux limiter used. See iMesser et al.l (|1998f ) and 
IPons et all (l2000h for discussions of the affect of different 
flux limiters. 

Models using the identical microphysics and initial con- 
ditions described in sect ion [5] were calculated using the 
EFLD code described in lRoberts et al.l (|2012f l with con- 
ve ction turned off. This code is similar to the one used 
in IPons et al. ( 1999) and is completely different from the 
one used in this paper. A derivation of EFLD in the con- 
text of this paper is included in Appendix [B] The EFLD 
luminosities as a function of time are shown in figure [H 
alongside the luminosities from section [5j EFLD clearly 
does a reasonably good job of predicting the total neu- 
trino luminosity, but poorly predicts the luminosities of 
each flavor. 

At early times, some deviation in the total luminosity 
is expected because the mantle, which is driving most of 
the neutrino emission, is not particularly optically thick. 
Additionally, at late times when the whole PNS becomes 
optically thin, EFLD deviates from the variable Edding- 
ton factor solution. But for the bulk of the PNS evolution 
the deviation between the two methods is around 10%, 
which is surprisingly good agreement. Given that most 
predictions made using EFLD codes have relied only on 
the total neutrino luminosity, it seems that previous re- 
sults can be reasonably trusted. The total luminosity 
emitted from the PNS is set at the neutrino spheres of 
each flavor, which is the last point at which EFLD can 
be considered reliable. The outermost layers of the PNS 
in which the neutrinos decouple can come into radiative 
equilibrium on a short timescale and therefore rapidly ad- 
just to the flux being pushed through them from below. 
As neutrinos propagate through the outer layers in the 
EFLD formalism, flux may be shifted between flavors un- 
realistically but the outer layers of the PNS will rapidly 
evolve to carry the right total luminosity. Therefore, it 
is not surprising that EFLD gets the total luminosity 
right but fails to predict the luminosities of specific fla- 
vors. Of course, EFLD makes no predictions regarding 
the spectral properties of the neutrinos. 

6.2. Neutrino Spectra and The Composition of The 
Neutrino Driven Wind 

The most striking difference bet ween the present sim- 
ulations and other recent studies (jHudepohl et aU\201Ck 



100 | i 




Time (s) 



Figure 6. Luminosities as a function of time for u e (solid black 
line), & e (dashed black line), v x (dot-dashed black line) of time 
and the total luminosity (red line) in the EFLD approximation. 
The gray and orange lines is the data from figure [2] The bottom 
panel shows the ratio of the total EFLD luminosity to the total 
luminosity calculated using the new code. 



IFischer ~et al. 2010) is the greater difference in the present 
study of the electron neutrino and anti-neutrino average 
energies at late times. There are a number of possible 
reasons for this difference. 

One is the initial model chosen for the PNS evolution. 
Rather than use an initial model from a sep a rate c alcu- 
lation of core - collap se, both IHudepohl et ali (|2010l ) and 
IFischer et ali (|2010l ) follow the entire evolution of the 
supernova. In so far as the initial models are similar, 
the two approaches should give the same answer. The 
initial model used here is somewhat dated and was cho- 
se n mainly to facilit ate the comparison with the work 
of IPons et al\ (|1999h . At early times the initial pro- 
genitor model will certainly affect the properties of the 
emi tted neutrinos sig nificantly, but after the first sec- 
ond [PcmsTeFaD (|1999h found that the evolution does not 
depend sensitively on the initial progenitor model. Of 
course, the difference in the average energies of the elec- 
tron and anti-electron neutrinos is a fairly subtle effect. 
Therefore, the effect of the progenitor model should not 
be ruled out, but based on the argument below it seems 
unlikely that the progenitor model is the dominant fac- 
tor. 

It is possible that the methods used for transport dif- 
fer enough to give disparate results. This also seems 
unlikely considering all three approaches come close to 
directly solving the Boltzmann equation, that the for- 
malism de scribed in this work is f airly s imilar to the for- 
malis m of IHudepohl et all (|2010D (see iRampp fc Jankal 
l2002f). and that the appro aches of IFischer et ali (|2010l ) 
and lHiidepohl et al\ (l2010f) have been shown to yield sim- 
ilar results (jLieb cndorfe r et al.1l2005f ). 

A more significant difference though may be the mi- 
crophysics employed. The difference between the elec- 
tron neutrino and anti-neutrino spectral temperatures 
is mainly set by the difference between their respective 



Proto-Neutron Star Evolution 



15 




Time (s) 

Figure 7. Energ y and number luminosities for a model using 
the IBru cnn ( 1985) nucleon capture rates compared to the fidu- 
cial model described in section[5] Notice the slightly increased v e 
and v £ cooling rates at late times and the convergence of all three 
lu minosities . The black lines and red lines are for the model us- 
ing Bruenn (1985) rates, while the gray and orange lines are for 
the model described in section \E\ In the bottom plot, the ratio 
JV„ e / Njj^ , which is of consequence to the electron fraction in the 
neutrino driven wind, is shown in the bottom plot of the second 
panel. 



mean free paths to capture on nucleons, as the scatter- 
ing mean free paths for both species are nearly equal. 
Due to de-leptonization, there are far more neutrons to 
capture electron neutrinos than protons to capture elec- 
tron anti-neutrinos. Of course, it is possible for both of 
these reactions to have strong final state blocking (elec- 
tron blocking for the neutrinos and neutron blocking for 
the anti- neutrinos). If it is assume d that there is no en - 
ergy transfer to the nucleons, as in Fischer et al.l (|2011f ). 
then both reactions will be st rongly blocked, the elastic 
interaction rates described in iBruennl (|T985) go to the 
same value, and it is expected that average electron neu- 
trino and anti-neutrino energies will be similar at late 
times due to the similar charged current mean free paths 
for both species. 

The final state blocki ng symmetry p redicted by the 
charged current rates of IBruennl ([19851 ) does not agree 
with more detailed calculations of the electron neu- 
trino capture rates. There is in fact a strong asymme- 
try between the two reactions, because there is signifi- 



cantly more energy available in the entrance channel for 
v e + n — > e~ +p than for D e +p — > e + + n. The difference 
between the energy of the entrance channels is just the 
difference between the fermi energies of the neutrons and 
protons. The Fermi energies for interacting nucleons are 
given by eF,i = k% i /2M{ + Ui, where Ui is an isospin 
dependent potential energy due to strong interactions in 
the medium. For neutron rich conditions, the neutron 
potential energy is larger than the proton potential en- 
ergy due to the nuclear symmetry energy. Most of the 
potential difference, Un — Up is transferred to the outgo- 
ing electron in the reaction v e + n — > e~ + p. This effect 
can significantly decrease the absorption mean free path 
for electron neutrinos. Due to the large value of the nu- 
clear symmetry energy relative to the value expected for 
free nucleons, Un — Up accounts for a significant fraction 
of fl. Although this amount of energy is often not enough 
to put the final state electron above the electron Fermi 
surface, it is enough to put the final state electron in a 
relatively less blocked portion of phase spac e. This effect 
is incl uded in the relativistic formalism of iReddv et al.l 
( 1998), which is used to calculate the neutrino interaction 
rates used in the models presented in this work. The de- 
tails of the importance of realistic kinematics on charged 
current rates will be discussed in future work. 

To illustrate how more realistic rates affect the pre- 
dicted neutrino properties, a model identical to the one 
described in section [5] was run, exce pt that t he nu cleon 
capture rates were replaced with the IBruennl (j 19851 ) cap- 
ture rates neglecting the nucleon potentials. The lumi- 
nosities as a function of time are shown in figure [7J The 
changes in the luminosity are relatively small. The most 
obvious difference is that the luminosities of all neutrino 
species asymptote to one another at late times , whic h 
is similar to the behavior seen in iFischer et al.l (|2011l ) . 
From one to ten seconds, there is a significantly smaller 
difference between the luminosities than in the fiducial 
model. Cooling via electron neutrinos and anti-neutrinos 
is also increased at late times, but the v x luminosity is 
virtually unchanged, as expected. Before 1 s, the elec- 
tron neutrino luminosity is reduced. It is unlikely that 
this is significant, as the first approximately hundred mil- 
liseconds of these simulations are suspect for the reasons 
described above. The reason for the early time decrease 
is less clear. It is unlikely that this is due to the effects 
of nuclear interactions because the region where electron 
neutrinos decouple is in the mantle which is at low den- 
sity. 

The evolution of the average neutrino energies are 
shown in figure [5] There is little change between the 
models in the electron anti-neutrino, fj,, and r neutrino 
average energies, but there is a significant change in 
the electron neutrino average energies. At early times 
the average energy is reduced compared to the fiducial 
model and at late times it is increased. The late-time 
convergence is easily explained by the argument given 
in the paragraphs a bove and by the arguments given in 
IFischer et all <[20TTh . 

This difference is important to the composition of the 
neutrino driven wind. The electron fra ction of the neu- 
trino driven wind can be estimated as (|Qian fc Wooslevl 
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Figure 8. Top panel: Energy moments of the outgoing neu- 
trino flux me asured at the surface of the calculation using the 
IBruennl 119851 ) approximation for the electron neutrino and anti- 
neutrino capture rates on nucleons (black lines). The gray lines 
are f or the fiducia l m odel using the full capture cross-sections 
from |Rcddy ct alj 11998T ). There is little variation for the 9 e 
and i/ x energies between the two cross-section prescriptions, but 
there is a significant change in the u e average energies. Bottom 
panel: Predicted neutrino driven wind electron fraction as a func- 
tion of time. The dotted line is from a PN S model using the 
IBruennl j!985ft rates, the solid line is for the [Reddv ct al] ^1998fl 
rates including tensor polarizati on corrections and the mean fields 
(Horowitz & Percz-Garci'a 2003), the dot-dashed line is a model 
using the Re ddv et alj 1199811 rates without tensor polarization cor- 
rect ions and with mean fields, and the dashed line is a model using 
the IReddv et al.l 119981) rates with tensor polarization corrections 
but neglecting the effects of mean fields. Note that neutron richness 
is predicted from about 1.5 to 10 seconds when realistic kinemat- 
ics is used in the capture rates, while the wind is predicted to be 
proton rich throughout when the effects of the neutron and proton 
potentials are ignored. 
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where (ct) are the energy averaged cross-sections for neu- 
trino capture on nucleons, which are approximately pro- 
portional to e 2 . Smaller relative v e average energies and 
lower de-leptonization rates lead to a lower electron frac- 
tion in the wind. 

The evolution of the electron fraction in the neutrino 
driven wind calculated using equation [M] for both mod- 
els, as well as a model that does not include weak mag- 
netism corrections and a model that does not include 
mean field effects, but which do include full kinemat- 
ics in the structure functions, is shown in the second 
panel of figure El The captu re rates for low densities 
given in iBurrows et al.1 (|2006l ). which include first order 



weak magnetism and recoil corrections, have been used. 
This was done to put the comparison between the models 
on even footing, although it is not necessarily consistent 
with the rates used inside the PNS itself. The alpha ef- 
fect (jFuller fc Meverlll995D has also not been taken into 
account, which will push Y e closer to a half in both pro- 
ton and neutron rich conditions. Energy moments of the 
neutrino flux are taken using the values at the surface of 
the computational domain, not at infinity. The N Vts /Np a 
term is increasing with time in both models (see figure[7]), 
which increases the electron fraction in the wind. Note 
that once neutrinos are free streaming, this term is in- 
variant with radius. 

With these assumptions, the fiducial model of section 
[5] actually does result in a period of ne utron richness in 
the w ind, i n contrast to the re sults of IHiidepohl et all 
(120101) and iFischer et al.l (|201 If ). The wind is not very 
neutron rich (Y e > 0.45 at all times) and this change, 
by itself, would not result in substantial r-process nucle- 
osynthesis in the stan dard neutrino drive n wind where 
entropies are < 150 (jRoberts et all 120101 ). If for some 
reason the entropy were higher though, the possibility 
of an r-process remains . In contrast, the model that 
uses the rate s of | Bruennl (|1985f ) and the model using the 
IReddv et abl (119981 ) rates without isospin dependent nu- 
clear potentials consistent with the underlying equation 
of state results in a wind that is always proton-rich. 

To emphasize that this result is mainly due to the re- 
action kinematics and not the inclusion of weak mag- 
netism, models with the tensor polarization set to zero 
are also shown in the bottom panel of figure |S1 As is ex- 
pected f rom the fi r st ord er weak magnetism corrections 
given in lHorowit9 ([2002D , allowing for a tensor portion 
of the response increase the difference between the elec- 
tron neutrino and anti-neutrino average energies. This 
results in a lower electron fraction in the case including 
the tensor polarization relative to the case without. Still, 
the change between these two models is only a fraction 
of the change in the electron fraction when the IBruennl 
( 1985J) rates are used. 

Given the sensitivity to the neutrino interaction rates, 
it is possible that further improvement of the treatment 
of electron neutrino and anti-neutrino capture will alter 
this conclusion in one direction or the other. Because 
the asymmetry between the electron neutrino capture 
rates depends on the value of the symme try energy and 
the sy mmetry energy varies with density ([Fattovev et all 
it may be that different nuclear equations of state 
alter the predicted neutron excess in the NDW. Varia- 
tions of the calculation of the rates, such as incorporat- 
ing the effect of correlations in the nuclear medium, can 
signifi cantly change the timescale of the neutrino emis- 
sion ([Reddv et aLlll999l ) and possibly the spectral prop- 
erties. Even if updated rates only change the rates at 
high density, this may affect N Ue / _/V Pc and thereby change 
the properties of the wind. Such extensions depend on 
the underlying equation of state, which is uncertain, and 
also on approximations inherent to many-body theories 
of strongly interacting systems. While this deserves fur- 
ther consideration, the results presented here seem to 
indicate that effects due to kinematics, degeneracy and 
mean fields are crucial. 

It b ears mentioning that the study of IHiidepohl et all 
( 2010) did allow for energy and momentum transfer be- 
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Figure 9. Spectrum of integrated neutrino emission over the du- 
ration of the simulation. Although the peak of the P e spectrum 
corresponds to the peak of the v x spectrum, the v x spectrum has 
a significantly harder tail. 



tween the nucleons and leptons ([Buras eFaT][2006l ). but 
did not account for the difference between the neutron 
and proton mean field potentials. Their rates were 
calculated within the ra ndom phase approximation of 
iBurrows fc Sawverl (|1999Q , which is an improvement over 
the mean field rates used in this work. Additionally, weak 
magnetism corrections were appro ximately incl u ded i n 
this study via the prescription of iBuras et al.l (2006). 
Their average energies were further apart than when the 
iBruennl (|1985l ) rates were used, but the difference was 
still not great enough to result in a wind with a neu- 
tron excess. The average energies of the electron neu- 
trinos and anti-n eutrinos also asymptot e to one another 
fairly rapidly in IHudepohl etatl (|2010n , in contrast to 
the present study. This is all reasonably consistent with 
the affect of neglecting the mean field potentials in the 
nucleon kinematics, but it is far from certain that this is 
the main source of discrepancy. Furthe r exploration of 
why t his work differs from the work of [Hudcpoh l et all 
(2010) is surely warranted. 

Even given these caveats, it is tantalizing that the wind 
is neutron rich in the fiducial model once again. Exten- 
sions of the standard neutrino driven wind model which 
include heating from a source besides neutrinos can pro- 
duce the r-process e ven for such modest neutron excesses 
(|Suzuki et al.l[2006T) . 

6.3. Time Integrated Spectra 

In figure |9l the integrated neutrino luminosity as a 
function of neutrino energy at infinity is shown for the 
model described in section [5] The time integrated av- 
erage energies of the neutrinos are (e Ve ) = 8.3 MeV, 
(e Pc ) = 12.2MeV, and {e u J = 11.1 MeV. Although the 
[i and t neutrinos are as hot or hotter than the electron 
anti-neutrinos at early times, the time integrated aver- 
age is weighted more strongly towards late times so that 
they in fact have a somewhat lower average energy. 

Time integrated neutrino spectra are in teresting for 
both nucleosynthesis via the ^-process (jHeger et al.l 
|2005| ) and for pre dictions of t he diffuse supernova neu- 
trino background (jAndol [20041 ). The neutrinos are non- 
thermal and are not easily described by an effective 
Fermi-Dirac distribution. Given the sensitivity of the 
^-process to the energy of the emitted neutrinos (espe- 
cially the energies above threshold) , it seems that calcu- 



lations of neutrino-induced nucleosynthesis needs to be 
done with more accurate neutrino spectra to check previ- 
ous results in the literature. These integrated spectra are 
only approximate however because a substantial fraction 
(20%) of the neutrinos are emitted during the first second 
of mantle contraction, this phase of evolution contributes 
the majority of the high-energy tail, and the mantle con- 
traction phase is most sensitive to the approximate initial 
conditions used. 

7. CONCLUSIONS 

A new code for following the evolution of PNSs has 
been described and some first results obtained. In section 
12.21 a formalism for moment based neutrino transfer with 
variable Ed dington f a ctors has b een described, ba sed on 
the work of lThornd (fl98l and iLindquistl (fl966l ). The 
framework is fully general relativistic and is formulated 
in the rest frame of the fluid, which simplifies calculation 
of the collision terms. The code employs energy inte- 
grated groups, rather than discrete energies, for solving 
the radiative transfer problem. This makes it well suited 
for dealing with problems were thermodynamic equilib- 
rium holds in large portions of the problem domain and 
the distribution functions may contain sharp Fermi sur- 
faces. A method for finding Eddington factors from a for- 
mal solution to the static Lindquist equation was also de- 
scribed. Additionally, general descriptions of the source 
terms for absorption, scattering, and pair annihilation 
have been provided which are consistent with the for- 
malism, explicitly obey detailed balance, and therefore 
naturally deal with the transition to equilibrium. The 
details of a fully implicit numerical implementation of 
these transport equations alongside the equations of gen- 
eral relativistic hydrodynamics were then described in 
section |4] 

The results of a fiducial model of PNS cooling were 
presented in section [SJ The evolution proceeds similarly 
to previous results in the literature in which a similar nu- 
clear equation of state and neutrino opacities were used. 
I have focused on the spectral properties of the emit- 
ted neutri nos, which were n ot well described by the for- 
malism of iPons et al. (1999). Similar behavior is found 
to other recent results in the literature: spectral soften- 
ing as a function of time and convergence of the v e and 
v TAL luminosities after about two second s of evolution 
(jFischer et aZ.ll20lot IHudepohl efaII[20Tl 

In contrast to other recent studies (c.f. iFischer et al.l 
1201 lh however, the new studies show that the average 
energy of the electron neutrinos does not converge to 
the average energies of the other neutrino flavors at late 
times. Additionally, the electron neutrinos are signifi- 
cantly cooler than the anti-electron neutrinos for most of 
the simulation. This difference is likely due to the treat- 
ment of charged current neutrino interactions, where a 
realistic treatment of the nucleon kinematics including 
the nuclear potential is important (see section 1572")) . The 
implications of this result for the electron fraction in the 
neutrino driven wind and possible r-process nucleosyn- 
thesis in this environment were discussed and warrant 
further exploration. 

A quantitative comparison was also made between the 
results of an EFLD calculation of PNS evolution and evo- 
lution with the code described in this paper in section lBTTl 
It was found that EFLD provides a good approximation 
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Figure 10. Comparison of H g / E g (crosses) and 31 (solid lines) for a purely absorbing (black lines) and an isotropic scattering sphere (red 
lines) with total optical depth one and a first order scattering sphere (blue lines) with t = 0.1. Aside from at the center, there is excellent 
agreement between the formal solution and the results of the moment calculation. The dashed lines are the second Eddington factors, 32 , 
for the same models. 



to the total neutrino luminosity during periods in which 
the neutrino luminosity is dominated by emission from 
optically thick regions. This approximation does break 
down in the optically thin regime as expected. Addition- 
ally, it does a poor job of predicting the luminosities of 
individual neutrino flavors. 

The most significant improvement which could be 
made to this work would be to improve the initial mod- 
els. This could be done either by using a more realistic 
post-core collapse initial models or updating the code to 
allow it to follow colla pse and bounce itself, sim ilar to 
iHiidepohl eta!\ (|2010D and iFischer et all (|2010D . Such 
improvements will affect the early time evolution, but 
are probably less important to the evolution of the PNS 
after one second. Additionally, a more realist ic equation 
of sta te that includes nuclei at low densities (jShen et al.1 
120111 ) should be employed with consistent opacities. 

Of course, this work has also been limited to one di- 
mension, which may be a gross, although necessary, over- 
simplification. Convection, magnetic fields, and rotation 



may be central players in the evolution of PNSs. Ap- 
proximate mixing length convection will be included in a 
subsequent version of the code. In the future, this code 
will be applied to understanding the diffuse supernova 
neutrino background, predicting the affects of different 
prescriptions for the nuclear equation of state on PNS 
cooling, and investigating black hole formation. 
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APPENDIX 
CODE TESTS 

Here I consider static transport through an homogeneous sphere with unit radius, and unit optical depth. A flat 
space-time is assumed. The first test performed is for consistency between the formal solution of the Boltzmann 
equation and the moment equations. In addition to the factors and 53, the formal Boltzmann solver can also 
calculate g\ = w 1 /w° which should be equal to H g /E g . For a purely absorptive atmosphere, the formal solution is 
exact. A comparison of g± and H g /E g is shown in figure [TO] For this calculation, one hundred equally spaced radially 
zones and a grid of 150 tangent rays with impact parameters spaced equally in radius were used. The calculation 
was then evolved for ten units of time. The differences between the two Eddington factors are negligible, aside from 
in the inner most zones. This agreement does not depend strongly on the number of tangent rays employed. The 
disagreement in the inner most regions is due to the small number of tangent rays which have impact parameters that 
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Figure 11. Top Panel: Eddington factors for radiation streaming from a homogenous sphere into free space. The solid lines are the 
Eddington factors gi for the purely absorptive sphere (black) and the isotropic scattering sphere (red). The crosses show H g /E g . The 
dashed lines are the second Eddington factor, 32, for the same models. The vertical dotted lines show the radii at which the opaque sphere 
ends. The horizontal dashed line shows the expected asymptotic value of 32 for free streaming radiation. Bottom Panel: Properties of 
the radiation field as a function of radius. The solid lines show the radiation energy density and the dashed lines show the luminosity per 
steradian, r 2 H g . The colors are the same as in the top panel. Once again, the vertical dotted line denotes the end of the opaque sphere. 



are less than the radius at which the Eddington factor is calculated. The distribution function is not well resolved and 
is therefore in error. Such problems should not arise in the actual evolution of PNSs, as the inner most regions are 
generally opaque. 

As a second test, a sphere which includes a scattering contribution to the opacity is considered. In this case, the 
formal Boltzmann solver no longer gives an exact solution of the transport equation because of the approximate 
treatment of the scattering terms (see section The problem set up involves a sphere of optical depth one with 10% 
of the opacity coming from absorption and 90% coming from isotropic scattering. Using a similar numerical set up to 
the absorbing atmosphere gives the results shown in figure 1101 Once again the deviation between the formal solution 
and the results of the moment calculation are small. Above the two innermost zones, the maximum deviation of the 
two results for g\ is less than 1%. At radii greater than 8, the deviation is less than 0.1%. A second scattering test 
problem was run with the isotropic scattering opacity, Xo> set to zero and the first order scattering opacity, xh se t 
to 90% of the total opacity. The sphere was assumed to have an optical depth r = 0.1. The agreement between the 
moment calculation and the formal solution was fou nd to be similar t o the i sotropic scattering case. 

Tests similar to the homogeneous sphere tests of iRampp fc Jankal ([2002) have also been performed. In these, a 
unit optical depth sphere of radius one is included inside a transparent region of radius ten, with a sharp transition 
region from the semi-opaque sphere to the surrounding vacuum. Such a scenario is similar to the neutrino decoupling 
region of PNSs, and is therefore an important test problem for any neutrino transport code for PNS evolution. The 
calculation domain is split up into 200 zones, equally spaced in radius and 301 a grid of 301 tangent is employed. Only 
the radiation density, E g and flux, H g , are evolved. The calculation is then run for fifty time steps, which allows the 
calculation to relax to steady state and forget the details of the initial conditions of the radiation field. Once again, 
one calculation was run with a purely absorbing opacity and a second was run with 10% absorbing opacity and 90% 
isotropic scattering opacity. 

The final Eddington factors and radiation field for these calculations are shown in figure QTJ For the purely absorbing 
calculation, the formal solution of the Boltzmann equation is exact. In the inner most zones, the distribution function 
is under resolved in angle due to the small number of tangent rays which pass through this region. This is not a 
problem for PNS simulations, as the optical depth in the interior is always much larger than one. Therefore, the 
Eddington factors are close to zero and have negligible affects on the moment transport solution. The deviation of the 
moment solution, outside the inner most region, from the formal solution is less than 1%. In the decoupling region the 
agreement is excellent. Additionally both g\ and gi asymptote to their expected values for free streaming radiation 
far from the core. 

EQUILIBRIUM FLUX LIMITED DIFFUSION 

In the diffusion approximation, the distribution function is approximated by f(u>, /i) = /o( w ) + /lO^)/- 4 - This results 
in 172 =93 — 0. The time dependence is then dropped and compression terms are then dropped from the first moment 
equation. Ignoring the terms containing derivatives with respect to energy (which will drop out in the end of our 
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analysis anyway), I find 



dw° 
dr 



4— w l 

or 



= 3s 1 



If the redshifted frequency is denoted as v — uje't', the redshifted source function is 
which gives the redshifted neutrino flux per energy 



wi = — T- 



dr 



4— w v 
or 



(Bl) 
(B2) 
(B3) 



where D(u>) = (I/A* + Xo ~ Xi/3) -1 - This amounts to assuming that the neutrino flux instantaneously equili- 
brates to the gradients in the neutrino number density. Because a number of terms have been dropped, this expres- 
sion does not guarantee that w i < wq, which means that it is possible for the neutrino fluxes to violate causality 
(|Levermore fe Pomr aning 1981). This problem is usually circumvented by introducing a flux limiter, which is a correc- 
tion to the dif fusion coefficient which depends on £ = w\/wo and serves to keep t he fluxes causal. In section the 
flux limiter of ILevermore fe Pomraning] (| 1 9 8 If) has been used, but it was found in lPons et al\ (j!999f ) that the results 
of EFLD calculations are reasonably insensitive to the choice of flux limiter. 
The total energy flux is given by 



H 



cLow 



-eT / duj 



dw° 
dr 



dr 



-w 



The equilibrium portion of the EFLD approximation constitutes assuming that fo(uj) 
equilibrium expression for w yields r r fit (,,\ 

The radial derivative is then given by 

d/c q (^) _ \ujdTe 



dr 



So that the energy flux is given by 



H 



T dr 
6tt 2 



D. 



dr 
dTe* 



/cqM(l ~ /eq(^)) 



(B4) 

f cq (uj,T, fx). Using the 
(B5) 

(B6) 



dr 



+ D 3 Te* 



dr) 
dr 



where the energy integrated diffusion coefficients are defined as 

/•oo 

D n = / dxx n D(xT)f cq (xT)(l - f cq (xT)). 
The number flux equation can easily be determined from the energy flux equation. This gives 



F 



6tt 2 



D 



dTef 
dr 



DoTe* 



dr) 
dr 



These expressions agree with the results of lPons et aH (|1999f K 

Now all that is left is to describe the evolution of the underlying medium. Using equation [ 
total internal energy of the medium including neutrinos is found to be 



dt n 



d 



da I £ 4 ^ 



2^ H 



(B7) 
(B8) 
(B9) 

the evolution of the 
(BIO) 



Using equation! 



the evolution of total lepton number is found to be 

§^(4^V^-^)=0. 



(Bll) 
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